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Section  (I) 
Introduction 


When  a  solid  body  moves  nonuniformly  through  a  compressible  gas,  compression 
and  expansion  waves  are  emitted  from  its  surface,  producing  a  reaction  on  the  body.  It  is 
of  interest  to  be  able  to  predict  the  pressure  forces  exerted  on  the  body  in  response  to  its 
motion.  In  general  this  problem  is  highly  nonlinear,  because  the  governing  equations  of 
a  compressible  gas  are  nonlinear.  If  the  interest  is  in  small  amplitude  perturbations  of  a 
rigid  body  in  uniform  motion,  then  the  governing  equations  may  rigorously  be  linearized 
about  the  steady  flow  produced  by  the  rigid  body.  The  deformation  induced  perturbations 
are  then  governed  by  a  linear  system  with  time  invariant  coefficients  (in  the  mean  body 
frame,)  even  though  the  mean  flow  field  itself  is  still  nonlinear.  Most  problems  of 
aeroelastic  interest  fall  in  this  category.  The  principal  significance  of  this  is  that  concepts 
of  generalized  force  derivatives  and  aerodynamic  transfer  functions  relating  generalized 
deformation  coordinates  to  their  generalized  forces  are  valid. 

The  problem  is  this:  there  are  no  good  methods  for  solving  the  linearized  equations 
of  motion  of  the  gas  significantly  more  efficiently  than  the  nonlinear  equations. 
Therefore,  there  is  little  advantage  (in  most  problems)  to  be  gained  by  linearizing. 

The  reason  for  this  is  that  the  linearized  equations  have  spatially  variable 
coefficients.  The  principle  of  superposition  applies,  but  there  are  no  known  solutions  to 
superimpose.  If,  however,  the  mean  flow  itself  is  uniform,  the  governing  equations  have 
constant  coefficients  and  elementary  solutions  are  known.  This  is  the  realm  of  classical 
linearized  aerodynamics,  or  equivalently  classical  acoustics,  for  which  quite  complicated 
solutions  can  be  constructed  simply  by  adding  together  many  simple  solutions. 

The  problem  of  constructing  time  dependent  flow  fields  governed  by  classical 
acoustics  will  be  reexamined  here.  It  should  be  kept  in  mind,  though,  that  the  underlying 
assumption...  namely  that  the  mean  flow  is  uniform...  is  fundamentally  flawed  unless  the 
mean  body  is  a  (generalized)  cylinder  moving  parallel  to  its  generators,  or  a  helical  sheet 
moving  tangent  to  itself,  these  being  the  only  rigid  bodies  that  can  slice  through  a  gas 
without  disturbing  it  Real  flight  vehicles  do  not  look  like  this,  and  do  produce  significant 
distortions  of  the  surrounding  gas.  For  slender  shapes,  typical  wings,  tails  and  fuselages, 
these  mean  flow  disturbances  ought  to  be  reasonably  small  except  near  stagnation  points, 
and  the  assumption  that  the  mean  flow  is  uniform  ought  to  be  reasonably  good  away  from 
stagnation  points,  at  least  in  terms  of  predicting  response  to  unsteady  deformations. 

Even  if  the  results  of  a  calculation  based  on  linear  theory  are  acceptably  accurate, 
though,  the  approximation  is  justifiable  only  if  those  results  can  be  obtained  considerably 
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more  quickly  than  they  could  be  from  a  full  nonlinear  analysis.  A  fast  approximate 
solution  can  be  useful  for  design  or  extended  parametric  studies  which  would  be  ' 
otherwise  prohibitive.  An  approximate  solution  that  is  not  faster  than  a  more  accurate 
method  is  useless.  It  is  not  clear  at  the  moment  whether  the  time  dependent  approach  to 
these  problems  is  in  the  "fast  useful"  or  "slow  useless"  category,  but  that  question  must 
eventually  be  answered.  For  the  case  of  rectilinear  motion,  the  equations  admit  simple 
harmonic  solutions,  so  that  the  problem  can  be  solved  for  a  single  frequency  of 
oscillation.  Methods  based  on  this  have  been  around  for  many  years  and  are  a  standard 
("fast  useful")  tool  in  aeroelasticity1 .  A  time  domain  method,  in  principle,  can  provide  a 
complete  spectrum  of  generalized  forces  from  one  transient  calculation,  or  be  coupled 
directly  to  the  structural  equations  of  motion  (including  nonlinearities)  to  look  at 
aeroelastic  response.  There  are,  then,  resons  to  believe  that  a  time  domain,  rather  than  a 
frequency  domain,  solution  method  may  be  useful. 

The  approach  that  will  be  adopted  in  this  report  is  to  start  from  basics,  the  inviscid 
equations  of  motion  of  a  gas  and  relevant  boundary  conditions.  Issues  of  linearization 
will  first  be  discussed.  Then  the  elementary  solutions  for  uniform  mean  flow  will  be 
described,  followed  by  a  presentation  of  the  general  method  of  superposition  of  those 
elementary  solutions  to  form  the  desired  solution  for  the  flow  about  a  given  deforming 
body. 

This  work  is  intended  as  an  extension  of  the  planar  wing  time  domain  analysis 
performed  earlier  by  M.  Blair  and  the  present  author2,5,  to  general  body  geometries.  The 
approach  taken  here  is  quite  similar  to  that  developed  by  Morino6,7,  but  with  some 
differences  in  the  integral  formulation. 
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Section  (II) 

Linearized  Inviscid  Equations 


The  governing  equations  of  fluid  mechanics  can,  for  our  purposes,  be  written  in  the 
inviscid  form  (Euler  equations): 


Dp/Dt  =  -pV*u 

n.i 

DuyDt  =  -— Vp 

n.2 

P 

Ds/Dt  =  0 

n.3 

where  D/Dt  is  the  material  derivative,  p  the  density,  p  the  pressure,  s(p,p)  the  entropy, 
and  ti  the  flow  velocity. 

The  relevant  boundary  condition  on  a  solid  surface  is  simply  that  the  fluid  shall  not 
penetrate  the  surface.  If  the  body  surface  velocity  at  some  point  is  a,,  the  unit  normal  to 
the  surface  at  that  point  is  n,  then  the  gas  velocity  at  the  surface  must  obey: 

u*n  =  u,*n  II.4 

These  equations  define  the  nonlinear  problem. 

Now  suppose  there  is  a  steady  flow,  denoted  by  subscript  0,  which  obeys  the  Euler 
equations,  and  which  is  perturbed  by  a  small  disturbance  flow,  which  we  denote  by 
subscript  1.  The  perturbation  flow  is  governed  by  the  linearized  Euler  equations: 


DotPi/pol =s— ~V*[pot?i] 

Po 

n.5 

D0U!  +u1-Vu0  =  — -i-Vp!  +  -^-Vpo 

Po  P5 

II.  6 

D0si  +ui’Vso  =0 

II.7 

where  Dq  denotes  the  linearized  material  derivative,  based  on  the  unperturbed  flow 

velocity  uo- 

These  equations  are,  of  course,  linear.  However  they  generally  have  variable 
coefficients,  which  depend  on  the  steady  flow,  and  are  not  significantly  easier  to  solve 
than  the  original  equations. 

If  the  mean  flow  is  uniform,  then  we  can  drop  all  terms  involving  gradients  of  the 
mean  flow.  Suppose,  in  addition  that  we  use  a  coordinate  system  which  is  fixed  relative 
to  the  undisturbed  gas,  so  that  the  mean  flow  velocity  is  zero.  Then  we  recover  the 
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classical  linear  acoustic  equations: 


Pi,t  =  -PoV-tfil 

n.8 

— *  1  T7 

Ui.t=— —  Vp! 

n.9 

Po 

o 

II 

*•* 

** 

mo 

Note  that  the  difference  between  classical  acoustics  and  linear  aerodynamics  lies 
not  in  the  governing  equations,  but  in  the  fact  that  in  typical  acoustic  problems  the 
sources  are  stationary  with  respect  to  the  gas,  while  in  aerodynamics  the  sources  move 
with  respect  to  the  gas. 

The  surface  boundary  condition  can  be  expressed  simply  as: 

no*ui=  IL11 

where  no  is  the  unit  normal  top  the  undeformed  body,  since  now  is  a  pure 
perturbation.  The  normal  velocity  of  the  body  must,  of  course,  be  small  compared  to  the 
speed  of  sound  for  the  linear  approximation  to  be  meaningful.  How  small  it  has  to  be 
depends  on  the  flight  Mach  number,  the  closer  the  Mach  number  is  to  1,  the  smaller  the 
normal  velocity  must  be  to  avoid  nonlinear  steepening  of  waves.  Note  also  that  there  will 
generally  be  some  points  on  the  body  where  the  local  body  velocity  is  along  the  normal. 
These  will  usually  be  local  points  of  failure  of  the  linear  theory. 

Now  it  is  apparent  that  in  this  approximation  the  entropy  perturbation  Sj  is  fixed  by 
the  initial  conditions  and  never  changes.  Hence  an  entropy  disturbance,  if  one  exists,  has 
no  effect  on  the  velocity  or  pressure  fields.  The  continuity  and  momentum  equations  can 
be  written  purely  in  terms  of  pressure  and  velocity: 

(Pi/Po),t  =  -a02V*u1  n.12 

ui.t  =  -V[p  i/p0]  n.13 

where  ao  is  the  undisturbed  speed  of  sound. 

We  can  further  split  the  problem  into  a  potential  wave  field  and  a  stationary 
rotational  flow  (gust).  To  do  this,  define  a  disturbance  velocity  potential  <p,  such  that: 

+  Pi/Po  =  0  11.14 

which  is  a  linear  Bernoulli  equation,  and  define  a  "rotational  velocity"  i?r  by, 

ur=u!— V«J>  11.15 

We  are  free  to  select  this  splitting  so  that  the  rotational  velocity  is  initially  solenoidal: 
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n.16 


V'u^.  =  0  at  t=0 

It  follows  immediately  that  the  rotational  part  of  the  velocity  is  independent  of 

time: 

^,t  =  0  H.17 

(and  is,  therefore,  always  solenoidal),  and  that  the  potential  obeys  the  classical  linear 
acoustic  wave  equation: 

□fy  =  V2<J>  -  <j>,t,t/a20  =  0  n.18 

Thus,  the  rotational  flow  is  determined  strictly  by  the  initial  conditions,  like  entropy. 
Unlike  entropy,  the  rotational  velocity  is  important,  though,  because  it  interacts  with  the 
wave  field  through  the  surface  boundary  condition: 

n<)*V<J>  =  n*a,-n0,u’r  11.19 

The  problem  that  must  be  solved,  then,  is  the  wave  equation  for  the  velocity 
potential  <J>,  subject  to  the  Neumann  boundary  condition  on  the  undeformed  body.  Once 
this  is  done,  the  pressure  can  be  determined  from  the  linearized  Bernoulli  equation.  Note 
though,  that  since  the  body  moves  with  velocity  T?s  the  time  derivative  in  Bernoulli’s 
equation  should  be  expressed  at  a  fixed  point  on  the  body: 

pi/po— (W-v*)  n.20 

where  <j>  is  the  time  derivative  of  surface  potential  at  the  point  on  the  body  which  moves 
with  velocity  "u,.  The  convective  derivative  can  be  decomposed  into  normal  and 
tangential  components,  and  the  normal  component  is  known  from  the  surface  boundary 
condition.  Hence,  it  is  sufficient  to  know  the  potential  on  the  surface  of  the  body  to 
determine  the  pressure  on  the  surface. 
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Section  (HI) 

Elementary  Solutions  of  the  Wave  Equation 


The  wave  equation  has  a  simple  spherically  symmetric  solution: 

0(x,t)  =  -Q(t)/(  4jiR  )  m.i 

where  R  =  x-xs  is  the  radius  vector  pointing  from  some  fixed  point  xs  to  the  field  point, 
R  is  its  length  and: 

t  =  t-R/ao  III.2 

is  the  time  at  which  a  sound  wave  must  leave  xs  to  arrive  at  x  at  time  t,  the  so  called 
"retarded  time."  This  solution  represents  the  potential  due  to  a  point  source,  emitting 
fluid  from  the  source  point  xs  at  the  rate  Q(t)  in  volume  per  unit  time. 

This  solution  is  valid  only  if  the  source  point  is  stationary,  that  is  if  xs  is 
independent  of  time.  If  the  source  moves  with  respect  to  the  gas  far  away  from  it,  the 
solution  is  modified  by  what  are  known  as  Doppler  effects: 

<(K^t)  =  -Q(T)/(47tRc)  m.3 

where  Rc  is  stretched  by  the  relative  motion: 

Rc  =  R|l-Mr|  m.4 

and  Mr  is  the  component  of  source  Mach  number  directed  at  the  receiving  point: 

M,.  =  R/R'M„  m.5 

Ms=^/ao  ni.6 

Note  that  if  the  source  Mach  number  component  is  1,  then  R«  is  zero  and  the  potential 
blows  up.  This  corresponds  to  the  field  point  being  on  the  Mach  cone  emanating  from  the 
source  when  the  source  moves  supersonically  through  the  gas. 

It  is  important  to  recognize  that  the  retarded  time  is  still  determined  by  the  same 
condition,  but  now  the  source  position  is  that  at  the  emission  time  x: 

t  =  T+  |x-x’s(x)|  m.7 

If  the  retarded  time  is  given,  the  reception  time  is  explicit.  If  the  reception  time  t  is 
given,  this  relation  is  a  nonlinear  algebraic  equation  for  x,  which  could  have  any  number 
of  roots.  If  the  source  Mach  number  is  always  less  than  1,  then  you  can  show  that  there  is 
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only  one  x  for  every  t 

The  source  solution  can  easily  be  used  to  generate  other  solutions  by  addition.  In 
particular,  since  any  derivative  of  a  solution  must  be  a  solution,  then  we  have  that: 

.<j>=V-[-K(x)/(47tRc)]  m.8 

which  is  the  potential  of  a  dipole  with  strength  and  orientation  fixed  by  K. 
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Section  (TV) 

Source  and  Doublet  Distributions 


We  first  consider  the  representation  of  an  arbitrary  wave  field  by  means  of  a 
distribution  of  sources  and  dipoles  over  a  given,  but  general  surface.  The  only  real 
difficulty  in  doing  this  is  that  we  must  adopt  some  method  of  parameterizing  the  surface 
over  which  the  singularities  are  to  be  distributed.  Since  the  retarded  time  depends  on  time 
and  position,  using  the  spatial  coordinates  of  the  surface  is  not  recommended.  Instead  let 
P0  denote  the  pair  (^,T|),  which  shall  uniquely  span  the  surface.  We  need  not  further 
specify  how  these  coordinates  are  to  be  defined,  except  that  they  shall  not  depend  on 
time.  Also,  let  dAo(Po)  denote  a  generalized  area  element  attached  at  the  point  Po-  This 
need  not  represent  a  physical  area.  In  fact,  we  may  simply  take  it  as  dAo  =  d^dq.  The 
presumption  is,  of  course,  that  the  coordinates  of  the  surface  are  specified  in  the  form: 

xJ)  =  xo  (P0,t)  IV.  1 

Having  adopted  such  a  parameterization,  it  is  then  apparent  that  the  potential  field: 

*(?.«=  -^JJ-3^^dA«(po) 

is  a  solution  of  the  wave  equation,  since  it  is  simply  the  sum  of  sources  q  dAo,  and 
dipoles  jxdAo.  The  reason  that  the  "area  element"  dAo  is  arbitrary  is  that  any  scale 
change  in  it  could  be  absorbed  into  a  rescaling  of  the  source  and  dipole  densities  q  and  |i. 
Note  that  at  each  integration  point  Po  the  retarded  time  must  be  found  from  the  delay 
equation: 

x  =  Hx-xofPo.'t)  |  /ao  IV. 4 

The  orientation  of  the  dipole  density  p.  is  arbitrary.  However,  if  it  were  tangential 

to  the  surface,  then  the  dipole  integral  could,  through  an  integration  by  parts,  be 

— > 

converted  into  a  source  distribution  and  absorbed  into  q.  It  is  sufficient,  then,  to  take  |i  as 
being  normal  to  the  surface.  Let  n  be  the  local  unit  normal,  so  that: 

5(P0,t)*n(Po,t)n(Po,t)  IV.5 

where  |i  is  the  magnitude  of  the  dipole  strength. 


rv.2 

IV. 3 
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Now  the  above  <j>  is  a  solution.  What  are  its  properties?  Consider  first  the  case  of  a 
pure  source  distribution,  ^u=0.  It  can  be  shown  that  die  potential  field  is  then  continuous 
but  that  the  normal  derivative  is  discontinuous  across  the  surface.  The  jump  in  normal 
derivative  determines  the  source  density: 

q(P0,t)  =  A(n*V<J))dA/dAo  IV.6 

where  dA  is  the  local  physical  area  elements,  and  A  is  the  jump  across  the  surface,  with 
positive  side  defined  by  the  direction  towards  which  n  points.  Note  that  since  the  jump  is 
determined  by  local  singularities,  the  retarded  time  does  not  come  into  play.  The 
problem  is  locally  incompressible. 

Now  consider  the  case  of  a  pure  dipole  distribution,  q=0.  It  can  be  shown  that  the 
normal  derivatives  are  continuous,  but  that  (since  the  dipole  term  is  the  derivative  of  a 
source  distribution)  the  value  of  potential  jumps  across  the  surface.  The  size  of  the  jump 
determines  the  dipole  strength: 

A4>dA/dAo  IV  .7 

It  is  evident,  then,  that  the  general  form,  containing  both  sources  and  dipoles,  will 
have  jumps  in  both  the  value  of  <j>  and  in  its  normal  derivative  across  the  surface,  and  that 
these  jumps  determine  the  source  and  dipole  densities  q  and  p.  If  the  surface  is  closed 
(and  its  exterior  is  of  physical  interest)  then  the  densities  q  and  p  have  no  direct  physical 
significance,  since  they  are  then  differences  between  exterior  and  interior  values.  If, 
however,  the  surface  is  open,  and  exposed  to  fluid  on  both  sides,  then  the  densities  are 
meaningful. 

Any  solution  of  the  wave  equation  for  the  wave  field  produced  by  a  body  can  be 
represented  by  some  distribution  of  sources  and  dipoles  over  its  surface,  as  described 
above.  But  while  arbitrary  choices  of  q  and  p  yield  a  solution,  that  solution  will  not 
generally  satisfy  the  flow  impermeability  condition  on  the  surface.  Our  objective,  then, 
is  to  determine  suitable  distributions  of  q  and  p  for  which  the  field  $  determined  by  them, 
satisfies  the  Neumann  boundary  condition  at  every  point  on  the  body.  There  are  generally 
many  ways  to  do  this. 

Note  that  an  open  surface  has  no  interior,  and  the  impermeability  condition  would 
usually  imply  equal  normal  velocities  on  either  side.  In  such  a  case,  the  source  density  q 
would  be  identically  zero,  and  p  would  measure  the  jump  in  potential  across  the  surface. 
This  is  the  model  of  a  zero  thickness  wing. 

For  a  closed  surface,  which  has  an  interior,  one  must  generally  make  a  choice  for 
either  q  or  p,  or  set  some  constraints  on  diem,  simply  because  there  is  one  boundary 
condition  at  each  point  on  the  surface,  but  two  degrees  of  freedom  (q  and  p). 
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Finally,  it  may  be  noted  that  x  the  foregoing  discussion  has  presumed  that  the 
sources  and  dipoles  are  to  be  placed  on^he  (undeformed)  surface  of  the  body.  This  is  not 
necessary.  In  fact,  singularities  can  clearly  be  placed  inside  a  closed  body  if  desired,  or 
the  entire  physical  surface  could  be  replaced  by  some  neighboring  simplified  simulacrum 
of  it  (in  fact  this  is  what  would  always  be  done,  the  only  question  is  to  what  degree  the 
model  is  faithful  to  the  reality.  Since  linearization  already  has  cost  dearly  in  fidelity, 
there  is  little  reason  to  be  too  religious  about  geometric  accuracy.) 
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Section  (V) 

Generalized  Green’s  Identity 


In  the  previous  section  we  looked  at  representations  in  terms  of  general  source  and 
doublet  distributions.  Usually  the  source  and  doublet  densities  involved  have  no  direct 
physical  significance.  In  this  section,  a  formal  derivation  of  Green’s  identity  for  an 
arbitrarily  moving  body  will  be  presented.  The  result  is  a  source/doublet  representation 
containing  only  the  known  normal  velocities  and  the  unknown  potentials  on  the  exterior 
surface  of  the  body. 

* 

We  begin  with  an  impulsive  source: 


G(x ,  t;  xq  ,  to  )  = 


8(t  —  to  +  R/ao) 
4tcR 


V.l 


where  5  is  the  Dirac  delta  function.  This  quantity  G  is  the  velocity  potential  of  a  "bang" 
source  emitted  from  =  "x  at  time  t  =  to.  It  is  a  Green’s  function  for  the  acoustic  wave 
equation,  in  that: 

tfG  .  V2G — GtUt  =  SOt-^o)S(t-to)  V.2 

a  o 


Note  that  since  the  source  is  impulsive,  the  field  is  zero  everywhere  except  on  the  surface 
of  the  sound  sphere  radiating  out  from  the  point  of  origin.  We  may  freely  think  of  the 
source  point  as  moving  around  in  space  along  any  path.  There  are  no  Doppler  effects 
because  the  sound  sphere  is  emitted  at  only  one  point  on  this  trajectory.  Any  motion  of 
the  point  when  it’s  not  emitting  can  have  no  influence  on  the  disturbance  generated. 

For  any  <|>  which  satisfies  the  wave  equation,  =  0,  and  any  G  that  satisfies  the 
inhomogeneous  wave  equation  Q^G  =  8,  the  following  is  a  simple  identity: 

V-(4>VG-GV<j))  =  <}>S+  -^[<(>G  t-G<f>it]>t  V.3 

ao 

Let  V  be  the  volume  exterior  to  some  closed  surface  S(t)  on  which  there  is  an  outward 
normal  n,  and  a  surface  velocity  u^(t).  Define  the  normal  component  of  the  surface 
velocity  as  Ujn  =  rr"us,  and  the  Mach  number  of  this  normal  velocity  as  =  u^/ao- 
Obviously  we  must  have  M®,  <1  or  the  linear  approximation  would  be  ridiculous.  We 
will  keep  terms  proportional  to  Mg,,,  though,  for  completeness. 

Integrate  the  above  identity  over  the  volume  V,  and  apply  Gauss’s  theorem  to  the 
divergence  and  Leibnitz’  rule  to  the  time  derivative  term.  We  get  the  following  identity: 
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$05>.t)S(t-lo)  -  J|dA[Gu„^-M»/a0(<|G1-G4,,)}-^  V.4 

A  =  (♦G.,-G*.tya02  V-5 

Un  =  n*V<j>  ;  Gn  =  n*VG  ;  V.6 

Here  Uq  is  the  normal  derivative  of  4>,  which  is,  of  course,  known  from  the  boundary 
condition,  u^  is  the  normal  component  of  surface  velocity,  which  is  also  known,  but  not 
the  same  as  Un  because  the  later  includes  the  effects  of  deformation  and  gust,  while  die 
former  includes  only  the  undeformed  body  motion. 

We  now  substitute  in  the  particular  choice  of  G,  and  group  terms  proportional  to 
5(t-to+R/ao)  and  its  derivative  8.  The  result  can  be  expressed  this.way: 

-4jc$(xo,t)8(t-to)  =  J|[8f0+8f1]+47t-^-  V.7 

where: 

v'8 

— > 

fl  -  Mm]  V.9 

1  aoR  R 

We  are  now  in  a  position  to  integrate  with  respect  to  time  t.  The  term  involving  A  will 
generally  produce  initial  condition  terms  in  the  solution.  These  are  of  no  interest,  and  we 
shall  take  this  contribution  to  vanish,  so  that  we  get  a  representation  of  $  at  the  arbitrary 
point  x$>,  tg  as: 

-4t«K4  to)  *  JdtJJ[dA[8fo+5fi]  V.10 

The  problem  is  that  S  depends  on  t,  so  we  cannot  simply  interchange  the  order  of 
integration.  Instead,  use  a  marker  parameterization  of  the  surface  Po  =  (£,il),  as  in  the  last 
section,  so  that  the  instantaneous  surface  S(t)  is  defined  by  the  collection  of  points 
xs(P0,t).  We  can  define  Jc(P0,t)  as  the  ratio  of  the  area  element  dA  to  the  generalized 
area  d^dt),  so  that: 

dA  =  Kd^dr)  V.ll 

Doing  this  allows  us  to  interchange  order  of  integration.  Now  we  must  perform  the  time 
integrals.  To  do  this  we  require  the  following  properties  of  the  Dirac  delta  function: 
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v  8(g  (t))=-^L  V.12 

lg(x)l 

S(g(t»= — §^E> —  V.13 

g(t)lg(t)| 

where  g(t)  is  any  regular  function,  and  x  is  a  root  of  g,  so  g(x)  =  0.  (If  there  are  multiple 
roots,  they  must  be  summed)  In  this  application,  we  have  g(t)  =  t-to+R/ao,  where  R 
depends  on  t  through  the  motion  of  the  surface  coordinates. 

Now  perform  the  time  integration  to  get: 

-4ic<|X*>,tb)  =  JJd^dn[Kfb-(0X+X4>)]/|g|  V.14 


$g 


The  integrand  is  evaluated  at  the  retarded  time  x,  the  root  of  g(t).  This  formula  is  still 
incomplete  because  it  contains  (in  fo)  the  partial  time  derivative  This  can  be  replaced 
with  the  total  time  derivative  of  $  at  the  body  fixed  point  Po  by  means  of  the  chain  rule 
identity: 

V.16 

where  we  have  decomposed  the  surface  velocity  T?s  into  its  normal  and  tangential 
components: 

Ui  =  uit+mijn  V.17 

The  end  result  for  the  potential  outside  the  body  can  now  be  written  in  the 
following  somewhat  more  compact  form: 

—4^00*0,  to)  =  JJ[CoUj,  +  Ci<j>  +  C2<j>  V.18 

Where  the  four  coefficients  Cj  depend  only  on  the  surface  geometry  and  state  of  motion, 
not  on  the  properties  of  the  wave  field  <|>.  These  coefficients  are  given  by: 
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„  Kd-M^2) 

00 — — 

Ci  =-<ra-f$/R2+X)/|g| 
r  _ 

3~  aoRc 

C2  =  -C3-V|g| 

X=-^-(n-^-Msa) 

aoRg 

g  =  l-u^*N/ao  =  1-H- 

R  =  xo-x(P0,x) 


Rc^RIgl 

N  =  R/R 


V.19 

V.20 

V.21 

V.22 

V.23 

V.24 

V.25 

V.26 

V.27 


The  above  formula  is  perfectly  general  (aside  from  having  ignored  initial 
disturbances  within  the  flow,)  and  makes  no  assumptions  about  the  shape  or  state  of 
motion  of  the  surface  S(t)  over  which  the  integration  is  performed.  It  is  one  possible 
choice  of  definition  of  the  general  souice/doublet  distribution  formulas,  which  has  the 
advantage  that  it  contains  only  physically  important  quantities,  the  normal  velocity  at  the 
surface  %  and  the  surface  potential  itself,  It  appears  more  complicated  than  the  earlier 
representation  only  because  the  derivatives  have  been  carried  through.  Some  small 
simplifications  are  possible  in  the  special  case  of  rectilinear  unaccelerated  motion  of  the 
surface  S,  but  they  are  insignificant  For  rectilinear  motion,  the  retarded  time  calculation 
can  be  performed  analytically.  For  other  types  of  motion,  it  would  have  to  be  done 
numerically,  though  the  additional  cost  would  likely  be  small  compared  to  the  total  cost 
of  solving  the  problem. 

Note  that  this  identity  applies  everywhere  in  the  exterior  space  surrounding  the 
body,  and  gives  us  the  potential  there  directly  if  the  normal  velocity  and  surface  potential 
were  both  known.  However,  only  the  normal  velocity  is  known  a  priori.  To  find  the 
surface  potential,  we  need  only  bring  the  field  point  xo  down  onto  the  surface  S(to),  at 
some  point  ^o,Tio»  say.  We  then  have  an  equation  mapping  the  surface  onto  itself: 

L*<p  =  b  V.28 

where  b  represents  that  part  of  the  integral  which  is  known  at  time  to,  and  L  is  a  linear 
Fredholm  operator  of  the  second  kind.  In  fact,  because  of  the  finite  time  delay  between 
points,  L  is  very  nearly  a  simple  scalar  multiplier. 
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Sec(V.l)  Retarded  Time 


The  retarded  time  x  is  the  solution  of  the  delay  equation  t=x+R/ao.  This  is  equivalent  to 
finding  the  positive  roots  of  the  equation: 

F(R)  =  R2- 1  xb-x(t-R/ao)  | 2  =  0  V.29 

Now  for  rectilinear  motion  this  is  just  a  quadratic.  Other  motions  would  be  at  least 
locally  rectilinear.  A  general  quadratic  iteration  is,  thus,  suggested,  based  on  the  Taylor 


series  expansion  about  some  point  R: 

F  =  F(R)  +  8R  F  R  +  .5[8R]2  F  r.r  V.30 

FR  =2[R-R-M,]  V.31 

F.R.R  =  2[1  -  M,2  +  R-M,/ao]  V.32 

Solving  the  quadratic  for  the  change  SR  gives  the  following  solution  for  R : 

R  =  [-B±V5]/A  V.33 

A  *  1  -•  M,2  +  R‘M,  /  ao  V.34 

B  =  R  (1  -  A)  -  R-M,  V.35 

D  =  (B  +  AR)2+A(|R|2-R2)  V.36 


The  two  roots  for  R  are  there  simply  because  the  equation  is  quadratic.  If  the  flow 
is  subsonic  there  can  be  but  one  positive  root,  and  the  other  is  simply  discarded.  At 
supersonic  speeds,  there  may  be  two  or  no  (real  positive)  roots,  depending  on  whether  the 
field  point  is  inside  the  domain  of  influence  of  the  source  point  or  not  If  there  are 
multiple  roots,  and  therefore  multiple  retarded  times,  the  contribution  of  each  to  the 
integral  is  added  in. 

In  rectilinear  motion  this  gives  the  exact  solution  regardless  of  the  initial  guess  for 
R.  In  nonrectilinear  motion  it  would  give  the  correct  result  if  iterated  starting  from  a 
nearly  correct  result  Since  any  real  calculation  would  involve  a  continuous  distribution 
of  retarded  time  over  the  surface,  one  would  always  have  a  good  initial  guess  in  hand 
from  the  last  point  processed. 
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Section  (VI) 
Uniform  Translation 


In  the  important  special  case  where  the  body  translates  with  constant  velocity, 
certain  simplifications  are  possible.  This  special  case  will  be  examined  herein,  first  for 
steady  flow,  then  for  unsteady  flow. 

Sec  (VI.  1)  Steady  Flow 


* 

Although  this  is  a  standard  problem  it  will  be  described  to  establish  the  method  and 
notation.  If  the  Mach  number  is  M,  directed  along  the  x  axis,  then  the  disturbance 
potential  obeys  the  classical  formula: 

<1-M2*,„-H|>,y,y-H|>.z.z=0  VL1 

If  we  define  p  -  VTl-M2),  and  do  a  transformation  of  variables: 

5p=x,y=Py,z=|3z  VI.2 

then,  as  is  well  known,  the  equation  reduces  to  Laplaces  equation: 

?Vo  VI. 3 

The  Green’s  function  is,  of  course,  -1/4x1^,  where  R«  is  the  radius  in  barred  coordinates. 
This  familiar  devise  is  the  steady  state  Lorentz  transformation,  which  we  will  make  use 
of  again  in  the  next  section. 

A  direct  application  of  Green’s  theorem  to  Laplace’s  equation  in  barred 

coordinates,  followed  by  a  transformation  back  to  physical  coordinates  in  the  integrals, 
yields  the  following  identity: 

— 4rc<J>(xo)  =  JJ[-j^-+P2  *  ]dA  VI-4 

Kc  Rc 

which  will  be  recognizable  as  a  distribution  of  sources  and  dipoles  over  the  surface. 
However  the  source  density  is: 

q  =  n*V<j>-M2<|>iXnx  VI.5 

which  is  not  known  from  the  flow  tangency  condition  on  the  solid  surface,  because  it 
contains  the  stream  wise  velocity  perturbation  $>x,  which  is  part  of  the  unknown.  A 
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similar  term  involving  tangential  velocity  appeared  in  the  general  unsteady  Green’s 
theorem  formulation. 

There  are  a  variety  of  ways  to  fix  this.  Here  we  use  integration  by  parts  to 
eliminate  the  tangential  velocity.  First  we  can  decompose  the  velocity  4>>x  into  normal 
and  tangential  parts,  with  the  result: 


<pPn2Un-M2  nx  yT*V<t>  VT.6 

A 

where  x  is  the  unit  tangent  vector  in  the  flow  plane,  and: 

Pn  =  Vi-Hi2  VI.7 

Mn  =  Mnx  VI.  8 

u„  =  n*V<j>  VI.9 

Y=i*t  VI.  10 

The  pan  of  q  which  contains  tangential  derivatives  can  be  integrated  by  pans  over 
the  surface.  We  suppose  that  the  surface  is  closed,  so  that: 

J/V-(M2nxY^/Rc)dA  =  0  VI.  11 

If  this  is  true,  then  we  obtain  the  representation: 

= //[pn2^-  +  <t>K  ]dA  VI.12 


where  K  is  a  somewhat  messy  function  of  position  and  Mach  number  given  by: 

K  =  K0/Rc3+K1/Rc  VI.  13 

Ki  =  (M2-3M2)/R,  VI.  14 

Ko  =  p2(l+Mn2)R-n  VI.  15 

The  variable  Rs  is  the  radius  of  curvature  of  the  surface  in  the  flow  plane.  It  appears  in 
the  result  because  of  the  differential  geometry  identity: 

V«(nxYt)  =  (l-3nx2)/Rs  VI.  16 

This  formula  contains  no  tangential  derivatives  of  the  unknown,  only  a  distribution  of 
known  sources,  and  a  weighted  integral  of  the  unknown  <j>.  It  naturally  contains  many 
terms  which,  if  they  were  not  negligible,  would  invalidate  the  assumptions  of  the  linear 
theory  on  which  the  formula  is  based.  (For  example,  you  must  have  Mn  «  1.)  This  is 
because  the  result  was  derived  for  an  arbitrary  body  shape,  while  not  every  body  shape  is 
really  suitable  for  a  small  disturbance  approximation. 
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Sec  (VL2)  Unsteady  Row 


Much  of  the  complexity  of  the  general  unsteady  formulation  arises  from  the  fact 
that  the  body  is  moving  relative  to  the  gas.  For  a  stationary  radiating  body  (acoustics), 
the  formulation  is  quite  simple,  and  gives  the  classical  result: 

-4**<&.t)  =  JJ[-JL^]dA+V0-/J[|^,t)]dA  VI.17 

where  Un  is  the  normal  velocity  at  the  surface  (known),  and  x  =  t-R/ao  is  the  retarded 
time. 

In  the  special  case  of  rectilinear  motion  of  the  surface,  we  can  transform  to  Lorentz 
coordinates,  which  make  the  wave  equation  identical  in  form  to  the  stationary  case: 

[^2--r-^rK>  =  0  VI.  19 

ag  3 F 

while  leaving  the  body  surface  stationary  (the  over  bars  denote  Lorentz  variables).  Now, 
since  the  bounding  surface  is  independent  of  time,  the  formal  solution  is  easily  obtained 
just  by  replacing  all  variables  in  the  stationary  body  integral  formula  with  Lorentz 
coordinates. 

The  result  is,  however,  awkward  to  deal  with  because  the  surface  coordinates  are 
deformed.  Therefore,  the  integrals  are  transformed  back  to  physical  space.  After  some 
tedious  manipulation,  we  get  an  integral  identity  with  striking  resemblance  to  the  steady 
state  formula: 

R  2 

-47t4x3))  =  u„  +  <J>  K  +  <j>  K2]dA  VI.20 

which  is  just  as  in  steady  flow  (K  is  the  same),  except  that  the  integrand  is  evaluated  at 
the  retarded  time,  and  there  is  a  new  term  directly  proportional  to  the  local  time 
derivative  of  <j>,  with  a  coefficient: 

a  _ 

,  nxY^Vt-»  „ 

K2  =  M2~ - R-n/Rc]  VI.21 

Note  that  if  the  Mach  number  is  set  to  zero,  the  formula  reduces  to  the  stationary 
radiating  body  case. 

This  formula  for  <J>  may  be  taken  as  the  basic  result  of  this  investigation.  It  give  the 
integral  representation  of  <{>  directly  in  terms  of  surface  values  of  <)>  itself  (with  no 
tangential  derivatives)  and  of  a  known  source  distribution.  It  is  little  harder  to  implement 
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than  a  steady  flow  solver  based  on  the  same  kind  of  approach.  In  fact  the  resemblance  is 
close  enough  that  the  best  way  to  implement  it  might  be  to  take  an  existing  steady  panel 
code  and  modify  it  to  account  for  the  retarded  time  and  the  extra  time  derivative  term  in 
the  equation. 


Sec  (VI.  3)  Discretization 

The  integral  representation  for  $  can  obviously  be  discretized  and  put  into  the  canonical 
form,  (on  the  surface  itself): 


[<M  =  Q[Un]r  +  Ci[<|>]r  VI.22 

where  the  brackets  [  ]  denote  vectors  spanning  the  surface  and  the  coefficients  Ck  are 
square  matrices  with  dimension  equal  to  the  number  of  elements  of  the.  surface 
representation.  The  values  of  these  coefficients  depends  on  the  particular  quadrature  rules 
adopted.  The  important  thing,  though,  is  that  for  the  case  of  rectilinear  motion,  these 
coefficient  matrices  are  independent  of  time.  In  particular,  Cq  and  Cj  are  the  same 
matrices  that  would  appear  in  a  steady  flow  solver.  (The  coefficient  matrix  C2  would,  of 
course,  not  appear  in  a  steady  solver  because  it  multiplies  a  term  which  vanishes  in  that 
case.)  They,  in  principle,  can  be  computed  once  and  stored.  The  only  drawback  to  doing 
this  is  that  for  a  large  problem,  these  storage  requirements  might  be  prohibitive.  If  there 
are  5000  elements  on  the  surface,  storage  of  all  three  coefficients  would  amount  to 
7.5  107  matrix  elements.  For  reasons  which  will  be  discussed  below,  only  a  small  part  of 
these  matrices  would  have  to  be  stored  to  do  a  direct  solution  for  [0]  by  time  marching. 

In  steady  flow,  the  C2  term  and  the  subscripts  "r"  can  be  dropped,  leading  to  a 
linear  algebra  problem  with  the  form: 

AIC[0]  =  Q>[un]  VI.23 

AIC  =  I-Ci  VI.24 

If  there  are  N  surface  elements,  then  AIC  is  a  full  N  by  N  matrix.  This  is  the  standard 
steady  flow  problem. 

In  an  unsteady  flow,  the  vector  [0]  indicates  the  current  time  vector  of  0’s  on  the 
surface.  The  vector  [0]r  (and  similarly  for  the  other  two  vectors  on  the  right  hand  side  of 
the  equation)  indicates  the  retarded  time  realization  of  0  at  each  surface  element.  The 
farther  away  the  sending  element  is  from  the  receiving  element,  the  farther  back  in  time 
the  signal  will  come  from.  Only  the  nearest  elements  would  depend  on  the  current  time 
(unknown)  solution,  so  that  the  size  of  the  linear  algebra  problem  to  be  solved  is  in  fact 
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far  smaller  than  the  number  of  elements  on  the  surface,  unless  the  time  step  is  so  large 
that  a  sound  wave  will  span  a  large  fraction  of  the  surface  in  one  step. 

To  make  this  concrete,  suppose  that  the  solution  for  $  is  stored  in  an  extended 
vector  [<j>]e  which  spans  the  surface  and  as  much  of  the  time  history  as  necessary  to 
include  all  possible  point  to  point  interactions.  (So  this  extended  vector  contains  copies 
of  [<|>]  at  equal  time  intervals.)  Then  by  the  retarded  vectors  we  must  mean: 

fcJr-HuJe  VL25 

[<Mr=TM>]e  V-26 

[<Mr  =  Ti[fle  VL27 

where  T  is  an  interpolation  operator,  and  Tj  is  a  difference  operator  (i.e.,  matrix 
representations  of  such  operators.)  For  example  we  could  use  Tx  =T(1-E)/At,  where  E 
is  a  unit  back  shift  and  At  is  the  time  step. 

If  there  are  N  total  surface  elements,  and  a  fixed  Nt  elements  within  the  discrete 
domain  of  dependence  of  any  given  element,  then  the  problem  can  be  reorganized  into  a 
simple  linear  algebra  problem: 

AIC  [<M  =  [b]  VI-28 

where  AIC  is  an  N  by  N  matrix,  with  at  most  Nt  nonzero  entries  on  each  row,  and  [b]  is 
an  N  vector  containing  all  of  the  known  data.  The  AIC  matrix  is  sparse,  but  its  sparsity 
pattern  is  unknown  and  would  vary  depending  on  how  the  surface  elements  were 
numbered  and  on  how  big  the  time  step  is  (and  so  on  how  big  the  domain  of  dependence 
of  any  element  is).  The  system  could  be  solved  by  a  one-time  LU  decomposition  of  AIC, 
followed  by  a  back  solve  at  each  time  step. 
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Sec  (VI.  4)  An  Exact  Solution 


One  difficulty  with  any  numerical  procedure  is  to  decide  whether  the  calculated 
results  are  correct  or  not.  For  the  disturbance  produced  by  a  general  three-dimensional 
time-dependent  deforming  body,  this  is  a  particular  problem,  since  exact  solutions  for 
typical  prescribed  body  motions  are  not  available.  This  difficulty,  however,  can  easily  be 
circumvented. 

In  the  Green’s  function  based  method,  the  surface  source  density,  which  is  simply 
the  normal  velocity  9<j>/5n,  is  set,  and  the  resulting  surface  potential  is  to  be  determined 
through  a  numerical  solution  of  the  integral  equation.  However,  consider  a  surface  with 
finite  volume  but  otherwise  arbitrary  shape.  Somewhere  in  the  volume  enclosed  by  this 
surface  place  a  moving  point  source  of  strength  c(t).  The  velocity  potential  induced  by 
this  source  is: 


<j>(x,  t)  =  -a(T)  /  4tcRc  VT.29 

From  this  elementary  potential  we  can  compute  the  velocity  induced  normal  to  the 
enclosing  surface: 

un  *  — ^-s-[  no*R(o/ao  +  o(  1  -M$  )/R<.  j-cno’Mo  ]  VI.30 

4x1^ 

^  — > 

where  no  is  the  unit  outward  normal  to  the  bounding  surface.  Mo  is  the  source  Mach 

number,  and  all  quantities  are  evaluated  at  the  retarded  time. 

If  we  set  the  normal  velocity  on  the  surface  according  to  this  formula,  then  the 

exact  solution  for  the  surface  potential  is  known:  it  is  just  the  simple  source  potential. 

The  advantage  of  this  ruse  is  that  it  is  applicable  to  bodies  with  completely  arbitrary 

shape,  moving  along  an  arbitrary  path.  Obviously  other  test  cases  can  be  built  similarly 

by  placing  more  sources  or  other  prescribed  singularities  within  the  body,  and  computing 

the  resulting  surface  normal  velocity.  The  simple  source,  however,  would  seem  to  be 

entirely  sufficient  for  checking  any  numerical  method. 

We  note,  for  future  reference  that  for  the  special  case  of  a  spherical  surface,  with 

the  source  placed  at  the  center  of  the  sphere,  the  surface  normal  velocity  simplifies 

slightly  to: 

un  =  [n0-Ro/a0  VI.31 
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Sec  (VI.5)  An  Example  Problem 


To  test  out  the  formulation  on  a  simple  example,  we  considered  the  following 
special  case.  Take  a  sphere.  Inside  the  sphere,  at  the  center,  place  a  simple  point  source. 
The  source,  and  the  sphere  surrounding  it,  translate  along  the  x  axis  at  some  fixed 
subsonic  Mach  number.  The  normal  velocity  on  the  surface  of  the  sphere  is  set,  as 
indicated  in  the  previous  section,  to  the  value: 

u„  =  [  no*R  o  /  ao  VT.32 

so  that  the  exact  solution  for  the  potential  should  be: 

<Kx,t)  =  -a(x)/4jcRc  VI.33 

for  any  arbitrary  choice  of  source  strength  a(t). 

The  problem  is  organized  in  the  following  way.  The  sphere  is  broken  into  N 
disjoint  panels.  The  distance  between  the  centers  of  panels  i  and  j,  Ry,  defines  a  signal 
time  delay: 

Ri,j  =  ao  At  (ny+ey  )  VI.34 

where  At  is  the  time  step,  ny  is  the  whole  number  of  time  steps,  and  £y  is  the  fractional 
part. 

The  surface  potential  is  discretized  as: 

-47t<t>i  =  Lj  [Co,yqjr]  VI.35 

where  the  superscript  "r"  quantities  are  defined  by  simple  linear  interpolation; 

qjr=q>Mq?^q?+1)  vi.36 

V  =  +eij  H>jli,+l )  VI.37 

<j>jr  =  (<j>^-(j>J«  +  1)/At  VI.38 

The  superscripts  njj  are  time  counters,  indicating  where  in  the  global  history  array,  the 
data  are  to  be  found  for  panels  i  and  j. 

At  the  start  of  each  time  step,  the  arrays  are  advanced,  and  the  current  time 
values  (J)1  are  set  to  zero.  The  potential  is  then  summed  and  the  result  of  the  summation  is 
stored  in  a  vector  bj. 

The  influence  coefficient  matrix  is  loaded  based  on  the  following: 
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VI.39 


AlCa  =  4jc 
if  ny  =  1,  then 
AIQj  =  AIQj  +  Q.yeij  +  Cz.y/At  VL40 

No  attempt  was  made  to  use  the  sparseness  of  AIC  because  the  problems  were  all  of 
small  size. 

The  influence  coefficients  were  evaluated  in  the  crudest  possible  way  by  a  simple 
midpoint  rule,  except  on  the  diagonals,  where  a  three-point  average  was  taken.  This 
method  is  unsatisfactory  and  was  used  only  to  test  out  the  algorithm. 

In  the  first  example,  the  Mach  number  is  0.5  and  the  source  strength  is  a(0=l-e_t, 
in  units  where  die  speed  of  sound  and  radius  of  the  sphere  are  both  1.  The  sphere  was 
divided  into  400  panels  by  20  equal  longitudinal  and  20  latitudinal  cuts.  The  time  step 
was  chosen  as  0.2.  The  result  for  the  time  history  of  the  potential  at  a  point  45°  over  the 
leeward  side  of  the  sphere  is  shown  in  Figure  1,  along  with  the  exact  solution.  The 
qualitative  features  of  the'  solution  are  reproduced,  but  the  numerical  results  are  rather 
jagged.  This  is  not  a  special  property  of  this  point,  as  shown  in  Figure  2,  which  shows  the 
longitudinal  surface  distribution  at  a  sequence  of  times.  It  is  believed  that  the  noise  in 
this  calculation  is  caused  primarily  by  the  poor  choice  of  near  field  influence  coefficients, 
not  by  the  low  order  temporal  interpolations  used.  The  fact  that  the  general  timing  of  the 
disturbances  is  correct  indicates  that  the  implementation  was  made  correctly.  All  that 
remains  to  be  done  is  to  replace  the  steady-state  influence  coefficients  by  more  accurate 
values. 

This  is  shown  somewhat  more  clearly  in  the  second  example,  which  is  identical  to 
the  first  except  for  a  smoother  source  distribution,  a(t)  =  (1-e-1)2.  The  distribution  of 
surface  normal  velocity  for  this  example  is  shown  in  Figure  3.  It  is  quite  smooth,  and 
approaches  a  steady  state  at  long  times.  The  corresponding  exact  solution  for  surface 
potential  is  shown  in  Figure  4,  and  the  numerically  computed  solution  in  Figure  5.  Notice 
that  the  numerical  solution  starts  out  correctly,  but  eventually  reaches  an  incorrect  steady 
state.  This  is  a  clear  indication  of  the  inadequacy  of  the  steady  influence  coefficients 
(which  are  what  determine  the  steady  solution.) 

The  numerical  part  of  this  work  was  not  carried  further  for  nontechnical  reasons. 
Much  remains  to  be  done,  principally  with  regard  to  the  implementation  of  accurate 
steady-state  influence  coefficients. 
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Section  (VII) 
A  Timing  Study 


An  effort  was  made  to  examine  the  cost  of  this  type  of  calculation.  Admittedly, 
such  costs  ate  very  machine  dependent,  but  the  results  are  suggestive.  The  major  part  of 
the  calculation  is  simply  the  accumulation  of  the  sums  representing  the  surface  integrals. 
Without  time  delays  this  amounts  to  simple  vector  multiplies  and  code  vectorization 
gives  a  factor  of  10  speed  up.  When  the  time  delays  are  included,  the  basic  operation  is  a 
summation  of  the  type  A(ij)X(k(i  j)),  which  does  not  vectorize  unless  the  vector  X  is 
gathered  and  scattered  before  and  after  the  operations  are  performed.  This  was  done,  with 
an  interesting  result.  The  calculation  became  a  trivial  part  of  the  calculation.  The 
gathering  and  scattering  operations  themselves  consumed  the  greatest  part  of  the  cpu 
time  (by  about  a  factor  of  8).  This  means  that  only  a  factor  of  2  savings  could  be  obtained 
by  vectorization,  overall,  since  the  gather/scatter  time  could  not  be  reduced. 

We  can  give  a  breakdown  of  the  costs  of  a  typical  time  step  for  each  of  six  distinct 
phases  of  the  vectorized  code.  The  units  are  percent  of  the  total  cpu  time  for  the  step: 

1)  evaluation  of  source  strengths  ....1% 


2)  source  terms .... 

3)0  terms  . 

....  gather  38% 

....  compute  10% 

4)  inversion  . 

. <1% 

5)  shift  stacks . 

. 1% 

It  is  apparent  that  80%  of  the  overall  cost  involves  no  calculations. 

It  was  estimated  that  in  vector  mode,  the  test  code  that  was  run  took  about  1.5 
10"3N2  seconds  of  CPU  time  per  time  step  (on  a  Gould  NP1  machine.)  Note  that  the  cost 
is  inherently  quadratic  in  the  number  of  panels  and,  therefore,  would  become  prohibitive 
rather  quickly.  With  this  run  time  estimate,  a  5000-panel  model,  run  through  1000  time 
steps  would  take  100  hours.  This  is  not  reasonable  given  the  physical  limitations  of  the 
underlying  theory.  The  Gould  is  about  1%  Cray2  speed,  depending  on  the  problem,  so 


24 


this  might  translate  into  1  hour  on  a  Cray.  It  is  quite  likely,  though  not  certain,  that  a 
direct  finite  difference  solution  to  the  same  problem,  with  the  same  resolution,  might  be 
faster  on  such  a  machine.  If  this  is  true,  then  methods  such  as  that  discussed  here  ought 
to  be  useful  only  when  N  is  small,  that  is  for  low  resolution  simulations.  Since  the 
resolution  of  the  underlying  theory  is  limited  anyway,  this  may  be  appropriate. 
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Section  (VII) 
Conclusions 


The  present  study  provides  a  general  formulation  of  the  time  dependent  linearized 
aerodynamic  problem  for  bodies  moving  along  arbitrary  paths,  and  in  particular,  for 
bodies  moving  rectilinearly  at  constant  velocity  (apart  from  small  amplitude  elastic 
deformations  of  the  surface.)  For  the  case  of  rectilinear  motion,  the  formulation  differs 
somewhat  from  previous  work  in  that  a  "direct"  Green’s  identity  is  given  which  contains 
only  the  surface  potential  and  surface  normal  velocity.  All  tangential  dipole  terms  were 
removed  by  an  integration  by  parts. 

A  significant  result  of  the  study  is  a  simple  exact  solution  that  can  be  used  to  check 
any  general  purpose  time  dependent  code,  at  least  for  bodies  with  finite  volume.  This 
simple  idea  could  be  useful  in  the  testing  of  other  codes,  though  it  does  not,  of  course, 
provide  results  for  problems  of  physical  interest. 

Implementation  of  the  proposed  formulation  in  a  working  code,  for  the  most  pan, 
remains  to  be  done.  One  example  calculation  was  performed,  using  rather  poor  near  field 
influence  coefficients,  to  test  out  the  basic  ideas.  The  results  of  that  calculation  indicated 
the  clear  need  for  a  better  near  field  representation,  if  accurate  solutions  are  to  be 
obtained,  but  gave  no  suggestion  that  there  would  be  any  difficulty  doing  so.  The 
example  did  allow  a  realistic  assessment  of  the  probable  computational  requirements  of  a 
full  simulation.  The  costs  would  be  quite  high  on  a  large  problem  (several  thousand 
panels),  though  of  course  no  method  is  likely  to  be  able  to  handle  a  complete  aircraft 
easily  and  cheaply. 
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a*  t/R 
Figure  1 

Computed  Potential  at  45  deg  leeward,  M  =  .5.o=l- 
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Figure  2 

Distribution  of  potential  with  longitude,  case  of  Figure  1 
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Figure  3 

Distribution  of  Normal  velocity,  M  =  .5,  a  =  (1  -  e_t)2 
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Figure  4 

Exact  distribution  of  surface  potential,  case  of  Figure  3 
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Figure  5 

Numerical  distribution  of  surface  potential,  case  of  Figures  3, 4 


